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Recently, the idea of classical Metropolis sampling through Markov chains has been generalized 
for quantum Hamiltonians. However, the underlying Markov chain of this algorithm is still classical 
in nature. Due to Szegedy's method, the Markov chains of classical Hamiltonians can achieve a 
quadratic quantum speedup in the eigenvalue gap of the corresponding transition matrix. A natural 
question to ask is whether Szegedy's quantum speedup is merely a consequence of employing classical 
Hamiltonians, where the eigenstates simply coincide with the computational basis, making cloning 
of the classical information possible. We solve this problem by introducing a quantum version of 
the method of Markov-chain quantization combined with the quantum simulated annealing (QSA) 
procedure, and describe explicitly a novel quantum Metropolis algorithm, which exhibits a quadratic 
quantum speedup in the eigenvalue gap of the corresponding Metropolis Markov chain for any 
quantum Hamiltonian. This result provides a complete generalization of the classical Metropolis 
method to the quantum domain. 
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Interacting many-body problems, classical or quan- 
tum mechanical, generally require an exponentially large 
amount of computing resources to find the (exact) so- 
lutions, as the system size increases. Nonetheless, in- 
genious classical methods such as Markov-chain Monte 
Carlo (MCMC) or quantum Monte Carlo (QMC) have 
been invented and proven to be highly successful in many 
applications. These methods, however, have certain lim- 
itations. For example, the running time of MCMC 
scales as 0(1/6) where 5 is the gap of the transi- 
tion matrix. For problems such as spin glasses where 
6 is small, MCMC becomes computationally inefficient. 
QMC methods, on the other hand, suffer from the neg- 
ative sign problem [2|. Despite many efforts have been 
made for improvement Q, this limitation is still one of 
the biggest challenges in QMC 

On the other hand, one of the most important goals in 
the field of quantum computation, as proposed by Feyn- 
man, is to look for new methods or algorithms that can 
solve these many-body problems more efficiently. This is 
referred to as quantum simulation, which can be imple- 
mented either by dedicated quantum simulators [U, or 
with universal quantum computers @. For the former 
case, high precision experimental techniques are required 
for faithful simulation, which are therefore closer to en- 
gineering problems. For the latter case, sophisticated 
quantum algorithms are needed. To this end, there are 
two main approaches, (a) bottom- up approaches (7- 12 1 : 
to look for entirely new algorithms based on the spe- 
cial properties of q uantum computers, and (b) top-down 
approaches T3-10|: to improve the existing classical al- 
gorithms by combining with elementary quantum algo- 
rithms. This work belongs to the latter class. 

Some years ago, Szegedy [3] described a method to 
quantize classical Markov chains. The key result is that 
a quadratic speedup 0(1/ VS) in the gap 6 of the transi- 



tion matrix is possible. This was later adapted to some 
new algorithms 16|, LL7[ that can prepare thermal states 
of classical Hamiltonians, based on the idea of quantum 
simulated annealing (QSA). 

In fact, Markov chains that correspond to classical 
Hamiltonians are relatively easy to construct, as all the 
eigenstates are in the computational basis |x). Prepa- 
ration of the thermal states of quantum Hamiltonians 
HI HI, IH, however, is a much more challenging prob- 
lem. It is because, classically, one needs to solve for the 
full eigenvalue problem (i.e., look for all eigenvalues and 
eigenvectors) for the quantum Hamiltonian first, which 
often requires more computational resources. 

The key to overcome this difficulty with a quantum 
computer is the use of quantum phase estimation al- 
gorithm the eigenvalues of a quantum Hamiltonian 
can be recorded without explicitly knowing the detailed 
structure of the eigenvectors. Based on this idea, Terhal 
and DiVincenzo [13( were able to extend the Metropolis 
algorithm to the quantum domain, avoiding the nega- 
tive sign problem in QMC. However, one limitation for 
their results is that the Metropolis step involve too many 
energy non-local transitions, making it conceivably inef- 
ficient. Recently, Temme et al., [H| addressed this prob- 
lem by introducing random local unitary operations in 
the update rule. In both cases, however, the underly- 
ing Markov chain is still classical in nature, which means 
that the scaling of the running time is still 0(1/ '5). 

In this letter, we propose a new quantum Metropolis 
algorithm for an arbitrary quantum Hamiltonian H at 
a given temperature T (j3 = l/fc^T), with a quadratic 
speedup 0(1/ V5). We call this algorithm the quantum- 
quantum Metropolis algorithm (Q 2 MA), as it shows a 
quantum speedup for the Markov chain of a quantum sys- 
tem being simulated. The previous version of the quan- 
tum metropolis algorithm [l|| is restricted by the no- 
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TABLE I: Comparison of various Markov-chain based algo- 
rithms for thermal state preparation 



Methods 

Q. Metropolis V 
Q. Metropolis IF* 
Q. Markov chain I 9 
Q. Markov chain II' 
Q 2 MA (This work) 



Hamiltonian Input b Output c Q. Speedup^ 



Quantum 


Any p 


pth 


No 


Quantum 


Any p 


pth 


No 


Classical 


| +) ®« 


CETS I 


Quadratic 


Classical 


| +) ®n 


CETS I 


Quadratic 


Quantum 


|a°) 


CETS II 


Quadratic 



"Quantum Hamiltonians include classical Hamiltonians. 

'Here, p is a density matrix, |+> = (|0> + |1)) /\/2, and |a°) is 
defined in Eq. J6]l. 

c Here p t ^ is the thermal density matrix (cf. Eq. Q), CETS II 
is defined in Eq. 1191 , and CETS I (see also Ref. [1JJ|) is similar 
to CETS II, but the {\i)} is replaced by the computational basis. 
Both of them are equivalent to p t ^. 

d We consider only the quantum speedup with respect to the gap 
8 of the transition matrix of the Markov chain. 

Tn Ref. [13H . the Metropolis rule is implemented by a controlled- 
swap with a "duplicated environment". 

■^In Ref. |lS(| . improvement of Ref. |l3|l , a rejection rule for quan- 
tum states is imposed to implement the Metropolis rule. 

9 In Ref. Il6ll . combining Szegedy's Markov-chain quantization 
with quantum simulated annealing (QSA). 

h In Ref. |l7l| . improvement of Ref. [Tr3l using Grover's fixed-point 
search. 



satisfies the detailed balance condition, 



(3) 



A solution which can satisfy the detailed balance condi- 
tion is rriji = SijZji, where Sij — Sji is any symmetrical 
transition probability, 



min{l, e 



} 



(4) 



is sometimes called the Metropolis filter. In the practi- 
cal implementation of the Metropolis method, one starts 
with some initial configuration, and then apply a random 
transition (e.g. single spin-flip) and compare the energy 
between the new eigenstate and the old one. If the new 
eigenstate has a lower energy, accept the move. Other- 
wise, accept the move only with a probability distribution 
given by the ratio of the corresponding Boltzmann fac- 
tors e~P( Ei ~ Ei \ This is called the Metropolis method. 

The performance of the Metropolis method depends on 
the properties of the underlying transition matrix M of 
the Markov chain, especially, the spectrum of the eigen- 
values Xk, which are all positive and bounded by 1, and 
the largest eigenvalue is always 1. For convenience, we 
order them as 



cloning theorem insofar that the information of a eigen- 
state cannot be retrieved after the proposed move in the 
Metropolis step. We relax this restriction by adopting 
a dual representation where the basis states consists of 
pairs of eigenstates related by time-reversal operation. 
For the cases where the physical system being simu- 
lated is time-reversal invariant, the pair of the eigenstates 
forming the basis vector become identical. A comparison 
of various Markov-chain based methods is summarized in 
Table H 

The goal of the Q 2 MA is to prepare the coherent en- 
coding of the thermal state (CETS) |a ) (cf. Eq. (JTSJl), 
which, after tracing out the ancilla qubits, is equivalent 
to 



Pth 



(1) 



where Z — Tr [e - ^] is the partition function, and \<pi) 
is the eigenstate of a quantum Hamiltonian H , associated 
with the eigenvalue E{. Before going into the technical 
details of this work, we briefly summarize the important 
features of Markov chains and Metropolis sampling, in a 
way applicable to both classical and quantum systems. 

Brief review of the Metropolis method — In the stan- 
dard Metropolis scheme, the Gibbs distribution e^ l3Ei /Z 
of certain eigenstates is generated through a Markov 
chain M, in which the matrix element my refers to the 
transition probability from the eigenstate \(pi) to \<fj). 
The equilibrium (stationary) distribution 



-pE t 



/Z 



(2) 



A = 1 > Ai > ... > A 



N-l 



> 



(5) 



The convergence time of a Markov chain is limited by the 
eigenvalue gap 5 = 1 — Ai of the transition matrix M by 
0{l/8) The purposes of this work are: (1) to show 
that the running time can be improved to 0(l/y/S), and 
(2) to fully extend the Metropolis sampling algorithm 
into the quantum domain. 

Generalization of the Markov-chain quantization - 
The ori gina l Markov-chain quantization method by 
Szegedy 14( is applicable to classical Hamiltonians only. 
To extend it to quantum Hamiltonians, we start with 
n qubits prepared in the state (l/-\/2) (|0) + \ l))® n , or 
equivalently N -172 ^^ \x), where N — 2". By per- 
forming a bit-by-bit CNOT gate on a set of n ancilla 
qubits initialized in the state |000...0), one can formally 
express the entire quantum state as 



JV-l 



(6) 



where \ipi) = J2x=o ( x I 1 * 9 *) \ x ) 1S an energy eigenstate of 
H (i.e., H \<pi) = Ei \ipi)), and 



N-l 



(7) 



is the time-reversal counterpart of \tfi), which is the 
eigenstate of the corresponding time-reversal Hamilto- 
nian H = H* with the same eigen-energy Ei (i.e., 
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H\(pi) — Ei \<Pi)). Through the phase estimation algo- 
rithm (PEA), the value of eigenvalue £7, can be obtained 
either from \tpi) or \tfi). This is the key property intro- 
duced in this paper which relaxes the constraints of the 
previous quantum Metropolis algorithm (l8j. In the fol- 
lowing, for the purpose of demonstration, we shall assume 
that the PEA can be applied perfectly, in the sense that 
each eigenstate can be uniquely identified by a unique 
eigenvalue. We leave our discussion on the effects of de- 
generacy on this algorithm in the Appendix. It turns out 
that the degeneracy of the eigen-energies of a quantum 
Hamiltonian alone is not necessarily a problem. 

Now, let us include into our Hilbert space an extra 
qubit initialized in |0) and define a more compact nota- 
tion, 



\i) = fa) m |o) 



(8) 



The information of a Markov chain can be encoded in a 
pair of unitary operators Ux and Uy (see Appendix for 
their detailed construction), 

Ux \i) = fa* \<P<) M |0) + 7,fc \<Pi) \v>k) |1» , (9) 
fc 

UY\j)=^2((Tjm\Vm)\Vj}\ty + 7jm \fjn) 1 1) ) , (10) 



where a 



Ik 



"kW 1 - z *k, and a kl 



(ipk \ K \ Here K is an unitary operator which plays 
the same role as the spin-flip in the classical Metropolis 
method, and Zik is the Metropolis filter defined in Eq. 
03} . Note that U Y is related to U x by a controlled-SWAP 
operation. 

Satisfying the detailed balance conditions — The de- 
tailed balance condition (Eq. can be recovered by 
the product of U x and U Y . From Eq. © and ([TO]), for 
j 7^ i, we have 

(j\U x U Y \i) = \a j f(z ij z ji ) 1 / 2 , (11) 

where we used (ifi\<Pj) = (<fj\<pi). On the other hand, 



(i\U x U Y \i)=\a u f+y2\a kl \ 2 (l- z ik ) 



(12) 



which, as we shall see, can be interpreted as the proba- 
bility of not undergoing a transition. 

Construction of the operator W — Now, using Eq. ([3]) , 
(fTT|) and ([H]), wc obtain the following decomposition: 

01 U X U Y \i) = 01 Dl'i \j) 0| M \i) (i| D-V* \i) , 

(13) 

where D n = Y)j—q \j) (j\ is a diagonal matrix, and 

M = m *j li) (*l> witn m « = (*l U X U Y \i) and = 
\aijj\ Zij for j ^ i is the transition matrix of the Markov 
chain. Within the subspace {\i}}, Eq. (Q~3|) implies that 
U x IIy and M are similar matrices, which means that 



they have the same set of eigenvalues (see Eq. (JSJ). 
Following Szegedy this property allows us to con- 
struct an operator 



W = (2A 2 - I) (2Ai - I) , 



(14) 



where 



JV-l 



Ai = 10 (*l and A 2 = U x U Y AiU Y Ux ■ (15) 

i=0 

The spectral properties of W can be seen in the following 
way: define \oth) = ^^q 1 a ki \i) to be the eigenvectors of 
K\U x Uy the eigenvalue equation can be written as 

AxU x U Y \a k ) = X k \a k ) . (16) 

On the other hand, using the fact that AiU x UyA-i — 
AiUyUxAi, we have, 



A 2 \a k ) = X k utu Y \a k 



(17) 



Eq. (fl6|) and (fT7| suggest that, if we start with 
vectors within Ai, W can be block-diagonalized into 
subspace of 2 x 2 matrices Wk spanned by the basis 
{\a k ) , U X U Y |a fc )}. Explicitly, 



cos(26» fc ) -sin(26» fe ) 
sin(26» fc ) cos(20 fe ) 



(18) 



where cos 9 k = Xk- Note that the eigenvalues of w k is 
e ±l0k . In the case of k — where Ao = 1 (or 8q = 
0), wo = I is simply an identity. From Eq. (| 1 3|) and 
the properties of the Markov matrix (see Appendix), the 
k = state is the coherent encoding of the thermal state 
(CETS) of Eq. p, 



\ a o) 



N-l 

E 

i=0 



(19) 



Recall that 



\ipi) |0), the state |ao) becomes 



the Gibbs thermal state pth = e~^ H /Tr{e~P H ] when the 
other qubits are traced out. 

One of the most important features about W is that 
the minimum eigenvalue gap A m j n = \29i\ of W is less 
than two times the square root of the gap S = 1 — Ai 
of the transition matrix M (using 29 > |l 
2Vl - cos 2 9): 



JliB I 



> 2V5 



(20) 



which is the origin of the quadratic speedup of Szegedy's 
algorithm. This completes our discussion on the neces- 
sary tools needed for the following discussion. 

Quantum simulated annealing ( QSA ) — Given a quan- 
tum Hamiltonian H and any finite temperature T, our 
goal is to obtain the corresponding coherent thermal state 
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of the form in Eq. ([19]) , from the initial state defined in 
Eq. (|6]), which can be readily prepared from the "all- 
zero" state |000...0), and be considered as the infinite- 
temperature state. To achieve this goal, we can use the 
method of quantum simulated annealing (QSA) . For 
completeness, we outline the basic strategy and summa- 
rize the related results in the following paragraph. 

The strategy of the QSA method is to prepare a se- 
quence, j — 0, 1, 2, .., d, of d+1 of coherent thermal states, 



N-l 



J\ - 



i=0 



(21) 



at a time; the temperature /3j = (j /d) ft of the coherent 
thermal state is lowered in each step. The basic idea is 
that, for sufficiently small A/3 = /3/d, one can show that 
\a J ) has a good overlap with ), i.e., 



IK + V >r >!- £ o 

where e is bounded by 0(A(3 2 \\H\\ 2 ) 
measurement 



AT-l 



n fc+1 = £ 



(22) 

then a projective 
(23) 



fe=0 



on the eigenbasis of 
temperature state \a 3 +1 ) 
whole process, 

I o\ ni I i\ n 

the total error (see Appendix) 

e = de < 0((3 2 \\H\\ 2 /d) 



to \a 3 ) will result in the lower 



with a high probability. For the 

o) , (24) 

(25) 



is suppressed by refining the step size, analogous to the 
quantum Zeno effect. 

Now, using the machinery we have developed, the 
operator Wj+i (the W operator defined in Eq. p^|) 
for |a^. +1 )) can be used to construct such a projective 
measurement, through the phase estimation algorithm 
[l7j (see also Appendix). Alternatively, one may per- 
form an artificial way of introducing decoherence |16| . 
where Wj+i is applied multiple times randomly. In any 
case, the number of controlled- Wj + i is at most 0(1/ V5), 
which is a quadratic speedup relative to classical Markov 
chains. This completes our description on the Q 2 MA. 

Conclusion — To summarize, we have described a new 
quantum Metropolis algorithm which extends Szegedy's 
method of classical Markov-chain quantization to the 
quantum domain, and provides a quadratic quantum 
speedup 0(1/ VS) in the gap S of the transition matrix 
M . The restriction encountered by the previous version 
of the quantum Metropolis algorithm is mostly due to 
the no-cloning theorem, where the required information, 



such as the associated eigenvalue, of an eigenstate cannot 
be retrieved after the proposed move in the Metropolis 
step. We relax this restriction by adopting a dual repre- 
sentation where the set of basis states consists of pairs of 
eigenstates related by the time-reversal operation. 

This result completes the generalization of the classical 
Metropolis method to the quantum domain. Morevoer, 
the advantages of this quantum algorithm over classical 
algorithms could be exponential, as there is no need to 
explicitly solve for the eigenvalues and eigenvectors in 
the classical ways for the quantum Hamiltonians being 
simulated. Finally, as the application of the Metropolis 
method to quantum Hamiltonians can be considered as 
a special case of quantum maps (operations) , it may be 
possible that the results presented here could be general- 
ized to allow quantum speedup for a much broader class 
of quantum maps. 
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Appendix: Construction of the unitary operators Ux 
and Uy 

Here we show how one may construct the unitary op- 
erator Ux defined in Eq. namely 

Ux\i)=^M<Pi)\Vk)\0)+lik\w)\<Pk)\l)) , (26) 

k 

where 

Cik = "fciV^fe > lik = C^kW 1 ~ Z ik ) (27) 

and 

a kl = {ip k \ K^) . (28) 

Here K is an unitary operator which plays the same role 
as the spin-flip in the classical Metropolis method, and 
Zik is the Metropolis filter defined in Eq. (j4j. Note that 
Uy is related to Ux by a controlled-SWAP. To ensure 
Eq. (|lll) is positive, we assume that K is symmetrical in 
the computational basis: 

(x'\K \x) = (x\K\x) . (29) 

For example, K can be the SWAP operation. 
We start with the following n-qubit state 

(i=y i (|0) + |l)f" , (30) 

which is equivalent to the "all-input" state 

JV-l 

-7= y>) , (3i) 

where N = 2". Suppose now we include a set of n ancilla 
qubits initialized in the "all-down" state 

|000...0) , (32) 

an apply a bit-by-bit CNOT operation, which is equiv- 
alent to a copy of the value of x to the register qubits 
(which is not the same as quantum state cloning), 

\x) |000...0) -> |a:) |a:) . (33) 

The resulting state is 

1 N ^ 

^ !>>!*> ■ (34) 



Given any Hamiltonian H, if we formally insert the 
completeness relation 

JV-l 

/=X><M¥>i| , (35) 

expanded in the eigenvector basis {|y>i)}, to the state in 
Eq. (f3"4"|) , then we get the state in Eq. © , 

1 w-1 

|«°>ee ^Mfa) - (36) 
without solving the eigenvalue equation, where 

JV-l 

ee £ falx) \x) (37) 

is the time-reversal counterpart of \(pi). This state can be 
considered as the infinite-temperature state, and is also 
the starting point for the quantum simulated annealing 
(QSA). 

We are now ready to consider the explicit procedure 
for constructing Ux defined Eq. ([9]). Starting with the 
paired state 

\<Pi) Wi) , (38) 
we apply the "kick" operator K to \<fii), and write 

K\<p i )=Y,<*#\<p k ) . (39) 

k 

Next, we implement the Metropolis filter by performing a 
controlled-rotation (based on the difference of the eigen- 
values): 

\ipi) \<p k ) |0) -> \ Vk ) (v^fc |0> + a/T - , (40) 
where 

Zij = {l,e-^- E ^} . (41) 

This creates a state as described by Eq. (|9]). Note that 
only the information about the difference of the eigen- 
values is needed. There is no need to determine each 
eigenvalue individually. 

Appendix: Generalization to include random "kicks" 

In applying the classical Metropolis method, for exam- 
ple to Ising model, one usual apply random spin-flips to 
the spins. This feature can be incorporated in our Q 2 MA. 
On the other hand, implementation of these spin-flips is 
also necessary for systems with time-reversal symmetry 
where the eigenstates \(pi) = \tpi) contain only real coef- 
ficients, in the computational basis. 
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To be specific, we consider a system of n spin-1/2 par- 
ticles, we include an extra ancilla qubits initialized as 
|000...0). Then, in the first step of Ux, we perform a 
transformation such that 



1000 



(42) 



Then, conditioned on each value of A, we apply a "kick" 
operator K\, where K\ = K\, e.g. spin- flip operator, to 
the spin-A. The net effect is that Eq. ^ becomes 

(j \U x U Y \t ) = -J2\(^\Kxm\ 2 V^ > ( 43 ) 



and Eq. (|10|) changes in a similar way. 



Appendix: Effects of degeneracy and the limitations 
of the algorithm 

The Q 2 MA presented in the main text does not neces- 
sarily break down when the Hamiltonian H of the quan- 
tum system is highly degenerate in the space of the eigen- 
states. For the sake of the argument, consider the Ising 
model: 



-ff ising = J ofcrj 

(id) 



(44) 



Although the eigenstates are highly degenerate, the 
Q 2 MA does not require the knowledge of the individual 
eigenvalues, but instead, it needs the difference between 
two eigenstates before and after the "kick" (spin- flip) . In 
the case of the Ising model, the change in energy is O(J), 
thus it is sufficient to ensure our resolution in resolving 
the energy change be smaller than O(J). 

To understand this point better, note that there are 
two places degeneracy would affect the agrument: (A) 
the implementation of the projector 



(45) 



in Eq. ([TS]). (B) The "leakage" to the degenerate sub- 
space in Eq. (J9j> . These two points are related. The pro- 
jector Ai can be implemented via a filter method (e.g. 
see Ref. [13]); this is essentially the same as applying 
the phase estimation algorithm (PEA) multiple times. 
Explicitly, we consider two states \<pi) \<fj), we wish to 



determine whether Ei 
w.r.t the operator 

U = t 



Ej. We have to perform PEA 



-iHt 



+iHt 



(46) 



which produces only the differences of two eigenvalues. 

Suppose we normalize all our eigen-energy such that 
< Ek < 1 .We determine an energy window A = 2~ a , 



for some integer a, the error e for a single run of the PEA 
is bounded by 



£ < 



A 2 



\Et-Ejl 



(47) 



If we apply PEA k times with the same precision, then, 

e -> e k , (48) 

which means that the errors for those energy changes 
\Ei — Ej\ being greater than A, |JSj — Ej\ > A, would 
become exponentially small. The remaining problem is 
to deal with those energies change smaller than A; this 
issue is related to the point (B) above. Let us call this 
new projector which has resolution up to A. 

Consider the transformation described in Eq. ([9]). 
Let us called the contribution from those Uik where 
\Ei — Ek\ < A a "leakage". It is because if those (nearly) 
degenerate eigenstates do not contribute to the resulting 
state, then the operator A[U X Uy^-'i would not trans- 
form any basis state \i) outside of the paired basis. 
In other words, ideally we want, for those cases where 
\E% — E k \ < A, 



(<fk\K \tfi) = (ideal). 



(49) 



Practically, this condition cannot be satisfied as typi- 
cally K should have a board or almost continuous 
spectrum. However, we can argue that the sum of these 
contributions, 



V ; 



£' 

k 



(50) 



could be made negligibly small for generic quantum 
systems where the Hamiltonian involves local interac- 
tion terms. Here the summation is over those k where 
\Ei-E k \ < A. 

The reason is as follows: let us focus again on Eq. ©. 
Before applying K, note that both \cpi) and \(pi) have 
the eigen-energy Ei. We may expect, in the worst case 
scenarios, \(pi) would contain a significant weight (i.e., 
sharply peaked) in those energy eigen-states close to \ 
In this case, the leakage in the degenerate subspace would 
be very bad. However, when we apply a "kick" K which 
does not preserve the symmetry of H , i.e., 

[K,H}^0 , (51) 

then the change in energy 

n^^KHK^-^H^i)] (52) 

is typically of order O(J), where J is the typical size 
of the local terms. If we set A <C J, then we expect 
that r\ <C 1, and the correction to Eq. (TTB1 and (|TT)) is 
0(r\). The random kicks described in the previous section 
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help spread out this effect (making the distribution more 
uniform). 

In short, as long as the distribution \ot k i\ = 
\{(fk \ K is smoothly distributed over a range of en- 

ergy which is much greater than the window A of the 
energy filter for A' 1; we should expect that the contribu- 
tion coming from the "leakage" can be made arbitrarily 
small, by decreasing A. We leave a more quantitative 
analysis of this point for the future work. 



improvement could be made by non-linear division of the 
Pj. Here we assume it bounded above, 



(^)<<^ 2 )„ • 



(57) 



For the whole process, the total error accumulates at each 
step to 

e = dxO(Ap 2 (H 2 ) )=O(p 2 (H 2 ) /d) , (58) 



Appendix: Error analysis on the process of quantum 
simulated annealing 

Here we perform an error estimation for the process 
of quantum simulated annealing (QSA). First, the pro- 
cedure of quantum simulated annealing starts with the 
infinite-temperature (/3 = 0) state Eq. ((6]), and end up 
at some finite-temperature (/3 ^ 0) state Eq. (fT9|) . The 
inverse temperature /3 is divided into uniform intervals 



(53) 



where j = 0, 1,2, ..,d, of d + 1. The coherent thermal 
states |oq) corresponding to the intermediate tempera- 
tures /3j are prepared sequentially. The is made possi- 
ble by the projective measurement which projects |q;q) 



to a; 



o 



The finally fidelity depends crucially on the 



overlap |(a^ +1 |a^)| 2 between these states. 
To estimate the overlap, note that 



as A/? — f3/d. Hence, the total error can be made arbi- 
trarily small by increasing d. In other words, to achieve 
any given accuracy e, one must perform at least 

d = O{p 2 {H 2 ) /e) (59) 
steps in the process of quantum simulated annealing. 

Next, for each step, the projective measurement 



n 



fc+i 



itf-1 

Ei 

k=0 



J+ 1 \ 



J+ 1 1 



(60) 



can be achieved by the phase estimation algorithm. For 
this purpose, an improved version is described in Ref. 17 1 
(Lemma 2). Here we summarize the result: to achieve 
an accuracy of eo i the number of application of the 
controlled- W gates is 



J+ 1 \ - 



-A/3H/2>i\ 



«n 



where A/3 = (3/d. The overlap, 



KK +1 )I 2 = 



-A/9ff/2|_i\|2 



("ole 



is second-order in A/3, i.e., 



J+l\|2 



1 -0(Af3 2 (H 2 )) 



(54) 



(55) 



(56) 



This is analogous to the the quantum Zeno effect. In 
general, the energy fluctuation (i? 2 ) is smaller for ther- 
mal states of lower temperatures. Therefore, a potential 



O (log (l/e ) /A min ) 



(61) 



where A m ; n is the minimun eigenvalue gap of W, see 
Eq. (|20|) . Here we should put cq — e/d, which equals 
e 2 //3 2 (H 2 ) fromEq. QJgfy. Recall that 1/A min is of order 
0(1/ VS) fromEq. (20]), we conclude that the number of 
controlled- W gate required is 



O 



(62) 



A quadratic speedup of 1/yS is achieved. This completes 
the error analysis. 



